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Abstract 

A  spectral  multidomain  penalty  method  model  has  been  developed  for  the  simula¬ 
tion  of  high  Reynolds  number  localized  stratified  turbulence.  This  is  the  first  time 
that  a  penalty  method,  with  a  particular  focus  on  subdomain  interface  treatment, 
has  been  used  with  a  multidomain  scheme  to  simulate  incompressible  flows.  The 
temporal  discretization  ensures  maximum  temporal  accuracy  by  combining  third  or¬ 
der  stiffly  stable  and  backward  differentiation  schemes  with  a  high-order  boundary 
condition  for  the  pressure.  In  the  non-periodic  vertical  direction,  a  spectral  mul¬ 
tidomain  discretization  is  used  and  its  stability  is  ensured  through  use  of  penalty 
techniques,  spectral  filtering  and  strong  adaptive  interfacial  averaging.  The  penalty 
method  is  implemented  in  different  formulations  for  both  the  explicit  non-linear 
term  advancement  and  the  implicit  treatment  of  the  viscous  terms.  The  multido¬ 
main  model  is  validated  by  comparing  results  of  simulations  of  the  mid-to-late  time 
momentumless  stratified  turbulent  wake  to  the  corresponding  laboratory  data  for 
a  towed  sphere.  The  model  replicates  correctly  the  characteristic  vorticity  and  in¬ 
ternal  wave  structure  of  the  stratified  wake  and  exhibits  robust  agreement  with 
experiments  in  terms  of  the  temporal  power  laws  in  the  evolution  of  mean  profile 
characteristic  velocity  and  lengthscales. 
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1  Introduction 


1.1  Background:  Localized  Stratified  Turbulence 


Stably  stratified  turbulent  flows  are  ubiquitous  in  natural  fluids  [15,32,52], 
Due  to  the  restoring  effect  of  the  buoyancy,  any  turbulence  producing  siz¬ 
able  vertical  transport  and  diapycnal  mixing  occurs  locally  and  episodically 
at  high  values  of  Reynolds  number  (defined  in  terms  of  a  characteristic  geo¬ 
metric  length  and  a  large-scale  velocity),  Re  >  O(104).  The  stable  stratifica¬ 
tion  gives  rise  to  characteristic  forms  of  fluid  motion  such  as  internal  gravity 
waves  and  remarkably  organized  clusters  of  quasi  two-dimensional  horizontal 
(“pancake”)  vortices  [16,44,49],  It  is  of  interest  to  fluid  mechanicists,  physical 
oceanographers  and  meteorologists  to  not  only  understand  the  fundamental 
physics  of  stratified  turbulence  but  also  correctly  parameterize  it  for  use  in 
large-scale  circulation  models  [36],  In  this  respect,  numerical  simulations  play 
a  critical  role. 

A  reliable  and  robust  numerical  simulation  of  stratified  turbulence  must: 

(1)  Employ  a  highly  accurate  discretization  scheme. 

(2)  Utilize  a  base  flow  which  retains  the  forcing,  boundary  conditions  and 
salient  kinematics/dynamics  of  geophysically  relevant  flows. 

(3)  Have  a  value  of  Re  as  close  as  possible  to  those  in  the  field  and  laboratory. 

(4)  Ensure  stable  long  time  integration  of  the  discretized  equations. 

It  is  well  known  that  spectral  discretizations  provide  the  most  accurate  nu¬ 
merical  approximation  of  the  governing  equations  [6,4,9],  particularly  in  DNS 
where  one  intends  to  capture  the  full  range  of  scales  of  the  flow  and  the  only 
dissipative  mechanism  at  the  smallest  scales  should  be  molecular  viscosity 
rather  than  some  truncation  error-induced  artificial  damping.  Stably  strat¬ 
ified  homogeneous  turbulence  is  an  extremely  popular  template  flow  for  a 
stratified  turbulent  flow,  in  part  due  to  the  ease  with  which  it  can  be  studied 
by  means  of  Direct  Numerical  Simulation  (DNS)  [30,35].  Periodic  boundary 
conditions  in  all  three  spatial  directions  may  be  considered  for  the  simulation 
of  such  a  flow  and  thus  one  may  employ  straightforward  Fourier  discretization. 
It  is  questionable  however,  as  to  how  representative  space-filling  stratified  ho¬ 
mogeneous  turbulence  is  of  its  localized  geophysical  counterpart. 


1.2  Representation/Discretization  of  Localized  Flows 


Introducing  a  non-periodic  vertical  direction  allows  for  the  specification  of 
a  more  realistic,  non-space-filling  stratified  turbulent  flow.  In  such  a  case, 
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Fourier  cosine/sine  discretization  may  be  used  in  the  vertical  to  maintain 
spectral  accuracy  [23,56].  Such  a  scheme  however,  is  restricted  to  strictly 
symmetric  Dirichlet/Dirichlet  or  Neumann/ Neumann  boundary  conditions  at 
the  top  and  bottom  of  the  computational  domain.  Low-order  finite  differ¬ 
ence  schemes  overcome  this  restriction  in  choice  of  boundary  conditions  but 
degrade  the  accuracy  of  the  numerical  approximation  due  to  their  inherent 
truncation  error.  A  logical  first  choice  is  to  then  use  orthogonal  polynomial 
base  functions  (Chebyshev  or  Legendre)  to  discretize  the  vertical  in  one  single 
domain,  as  has  been  done  in  the  simulation  of  a  stratified  shear  layer  [7] .  This 
approach,  however,  is  associated  with  a  fairly  slow  computation  when  any  siz¬ 
able  number  of  vertical  grid  points  (>  O(102))  is  desired. 

Spectral  multidomain/spectral  element  discretizations  [22,9]  allow  one  to  over¬ 
come  many  of  the  limitations  discussed  above.  The  computational  domain  is 
partitioned  in  subdomains  of  varying  size  and  order  of  polynomial  approxi¬ 
mation.  Subdomains  communicate  with  their  neighbors  via  a  simple  patching 
condition  [4],  The  immediate  advantages  of  such  an  approach  are  [28]: 

(1)  Flexibility  in  local  resolution  which  allows  treatment  of  strongly  localized 
solutions/energetic  flow  regions  without  overresolving  smooth/less-active 
parts  of  the  flow. 

(2)  Reduced  operation  count  and  impact  of  round-off  error  compared  to  a 
single  domain  computation  with  the  same  number  of  grid  points. 

Differences  in  the  formulation  between  the  spectral  multidomain  and  spectral 
element  methods  are  discussed  in  various  textbooks  [6,4], 


1.3  Spectral  Element  Simulations  of  High  Re  Flows 


Spatially  continuous  spectral  element  methods  have  been  applied  extensively 
to  DNS  of  engineering  flows  in  complex  geometries  [9],  which  becomes  costly 
when  even  moderate  Re  are  desired.  However,  as  mentioned  above,  geophys¬ 
ical  turbulence  tends  to  be  governed  by  high  values  of  Re.  The  efficient 
realistic-cost  simulation  of  a  high  Re  flow  is  inherently  under-resolved  posing 
significant  challenges  to  the  stability  of  a  spectral  element  (or  multidomain) 
scheme.  Addressing  these  challenges  is  an  ongoing  effort.  A  spectral  element 
ocean  model  has  been  used  in  the  the  simulation  of  either  the  global  oceanic 
circulation  [33]  or  coastal  processes  [34]  by  solving  the  two-dimensional  shal¬ 
low  water  equations  or  the  three-dimensional  hydrostatic  primitive  equations, 
respectively.  Under-resolution-driven  instability  is  overcome  in  part  through 
spectral  filtering,  the  application  of  an  explicit  low-pass  filter  on  the  spectral 
approximation  of  the  numerical  solution  [41],  However,  full  numerical  stabil¬ 
ity  is  only  ensured  through  introduction  of  enhanced  numerical  viscosities. 
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When  enforcing  the  no-slip  condition  at  the  domain  boundaries,  the  numer¬ 
ical  viscosities  generate  thick  and  non-physical  boundary  layers  which  may 
spuriously  bias  the  internal  high  Re  of  the  flow  [5].  In  the  case  of  the  in¬ 
compressible  Navier-Stokes  equations,  spectral  multidomain  techniques  have 
only  once  been  applied  to  localized  stratified  turbulence,  where  4th/5th  order 
subdomains  were  used  [54],  At  the  higher  Re  ~  O(104)  considered,  high  res¬ 
olution  (2562  x  512)  and  an  overlapping  subdomain  scheme  in  the  energetic 
core  of  the  flow  were  adequate  in  dealing  with  uncler-resolution-generated  in¬ 
stabilities  at  the  subdomain  interfaces.  The  only  other  existing  application  of 
higher  order  methods  to  moderately  high  Re  stratified  turbulence  have  been 
been  4th  order  compact  finite  difference  schemes  in  a  single  domain  [46].  Un¬ 
stratified  higher  Re  studies  using  spectral  elements  have  been  performed  using 
either  spectral  filtering  [43,17]  or  the  spectral  vanishing  viscosity  formulation 
[37,39].  Again,  the  Re  considered  are  relatively  moderate  (e.g.  cylinder  wake 
Re  =  3  x  103  or  channel  flow  of  Re  =  7.5  x  103).  The  spectral  element  method 
has  also  been  used  in  conjuction  with  the  dynamic  estimation  model  [1]  for  the 
simulation  of  Rer  <  103  channel  flow.  The  eddy  viscosity  coefficients  of  the 
dynamic  subgrid  scale  (SGS)  model  suppress  any  under-resolution  problems 
by  treating  the  energy  flux  to  the  unresolved  scales  as  a  dissipative  process. 
The  combination  of  such  dissipative  models  with  low-order  finite  difference 
schemes  in  the  non-periodic  vertical  direction  is  responsible  for  the  stable 
Large  Eddy  Simulation  (LES  )  simulation  of  higher  Re  of  channel  flows  [20] 
and  wakes  [12],  Although  the  truncation  error  of  these  low-order  schemes  pro¬ 
vides  additional  artificial  dissipation  and  thereby  achieves  higher  Re ,  it  also 
is  detrimental  to  the  accuracy  of  the  results. 

Following  the  above  discussion,  a  spectral  multidomain  scheme,  despite  its 
potential,  is  not  foolproof  to  numerical  instability  when  directly  applied  to 
a  strongly  under-resolved  simulation  of  an  incompressible  high  Re  flow.  The 
absence  of  any  artificial  dissipation  in  the  scheme  leads  to  instabilities  due  to 
aliasing  effects  [22]  closely  linked  to  under-resolved  numerical/physical  bound¬ 
ary  layers  [9].  Specifically,  Gibbs  oscillations  develop,  which  are  most  pro¬ 
nounced  at  the  physical  boundaries  and  subdomain  interfaces.  The  oscilla¬ 
tions  are  worse  with  increasing  Re  (and  thus  degree  of  under-resolution)  and 
generate  catastrophically  interacting  artificial  internal  waves  in  a  stratified 
flow. 


1.4  Spectral  Multidomain  Penalty  Methods 


Spectral  filtering  [41,22]  is  least  effective  at  the  boundaries  or  subdomain  in¬ 
terfaces.  Thus,  a  methodology  is  needed  to  ensure  the  stability  of  a  spectral 
multidomain  scheme  in  these  regions  for  a  highly  under- resolved  simulation.  In 
this  respect,  spectral  penalty  methods  are  very  effective.  Combined  with  spec- 
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tral  filtering  they  allow  the  highest  attainable  Re.  Penalty  methods  [27-29] 
recognize  that  numerical  instabilities  arise  because  boundary/patching  con¬ 
ditions  are  explicitly  enforced  and  there  is  no  provision  that  the  equation 
solved  is  satisfied  arbitrarily  close  to  the  boundary/subdomain  interface.  By 
collocating  a  linear  combination  of  the  equation  and  boundary/patching  con¬ 
ditions  (the  latter  multiplied  by  a  penalty  term)  at  the  corresponding  spatial 
locations,  the  penalty  method  produces  a  smooth  numerical  solution  with 
near- negligible  error  at  the  boundaries  and  interfaces.  Unlike  conventional 
spectral  element /multidomain  schemes,  penalty  methods  are  discontinuous  at 
the  subdomain  interfaces.  Spectral  multidomain  penalty  methods  (hereafter 
referred  to  as  SMPM)  using  moderate  spectral  filtering,  have  attained  a  Re 
value  as  high  as  3  x  105  [13].  Penalty  methods  have  been  used  successfully 
in  the  multidomain  simulation  of  turbulent  reacting  flows  [28,13]  and,  in  a 
geophysical  context,  of  the  shallow  water  equations  [21],  However,  all  the 
above  studies  deal  with  hyperbolic  partial  differential  equations.  Only  simple, 
time-independent,  one- dimensional  elliptical  problems  were  considered  in  the 
early  work  of  Funaro  [18,19].  The  only  work  which  has  implemented  penalty 
methods  in  the  simulation  of  incompressible  Navier-Stokes  equations  (specifi¬ 
cally  in  the  velocity-vorticity  formulation)  [53]  utilizes  the  spatially  continuous 
spectral  element  discretization,  employs  penalty  terms  only  at  the  physical 
boundaries  and  is  restricted  to  a  fairly  low  Re  =  103. 


1.5  A  New  Application  of  SMPM  to  High  Re  Turbulent  Flows 


This  paper  outlines  the  main  components  of  a  spectral  multidomain  penalty 
method  developed  for  the  purpose  of  highly  accurate  simulations  of  high  Re 
localized  stratified  turbulence  with  a  non-periodic  vertical  direction.  This  is 
the  first  time  multidomain  penalty  methods  have  been  applied  to  the  simula¬ 
tion  of  incompressible  flows,  with  a  particular  focus  on  subdomain  interfaces. 
In  the  numerical  method,  the  penalty  technique  is  combined  with  moderate 
spectral  filtering,  adaptive  subdomain  interfacial  averaging  and  a  high-order 
temporal  discretization  [38].  The  stratified  turbulent  wake  of  a  towed  sphere  is 
used  as  base  flow  for  model  validation.  The  SMPM  solver  replicates  extremely 
well  the  laboratory  wake  dynamics  both  in  terms  of  quantitative  timeseries 
of  mean  flow  quantities  as  well  as  qualitative  features  of  the  vorticity  held. 
The  transition  points  between  different  regimes  in  the  evolution  of  a  stratified 
wake  [47,49]  are  reproduced  correctly. 

The  basic  incompressible  stratified  how  model,  base  how  and  initialization 
procedure  are  discussed  in  §2.  The  numerical  method  is  presented  in  detail  in 
§3.  Results  validating  the  spectral  multidomain  method  are  shown  in  §4.  A 
summary  is  given  in  §5. 
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2  Incompressible  Stratified  Flow  Model 


2.1  Governing  Equations 


This  study  considers  incompressible  stratified  flow  governed  by  the  Navier- 
Stokes  equations  under  the  Boussinesq  approximation  [40,56]: 


<9u 

~dt 


X  1 

—  -[u  ■  Vu  +  V(u  •  u)]  +  F g - V/b  +  v  V2u 

1  _ ✓  Po  77b 

L(u) 


NW 


dp' 

~dt 


-V-(u(p’  +  p(z)))  +  kV2P'  , 


V  ■  u  =  0. 


where 


F„  =  -n—k 
P  o 


(1) 

(2) 

(3) 

(4) 


The  five  unknowns  to  solve  for  are  the  velocity  vector  u  =  ( u,v,w ),  and 
the  pressure  and  density  perturbations  p'  and  p',  respectively.  The  non-linear 
term  in  the  momentum  equations  (1)  is  written  in  the  skew-symmetric  form 
to  minimize  aliasing  effects  in  the  numerical  solution  [9].  The  perturbations 
p 1  and  p'  originate  from  the  decomposition  of  the  corresponding  total  values 
into  [40]: 


p  =  p(x,y,z)  +p\x,y,z,t)  ,  (5) 

p  =  po  +  p(z)  +  p'(x,y,z,t)  .  (6) 

Following  the  Boussinesq  approximation,  the  reference  pressure,  p(x,  y,  z)  and 
density,  p0  +  p(z)  are  in  hydrostatic  balance: 


dp 

dz 


“(Po  +  P)  9 


(7) 
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2.2  Boundary  Conditions 


The  computational  domain  is  designed  to  simulate  a  portion  of  a  laboratory 
water  tank.  The  top  and  bottom  boundaries  of  the  domain  correspond  to 
the  top  free  surface  (considered  to  be  free  slip)  and  bottom  wall,  respectively 
of  the  tank.  The  lateral  domain  boundaries  are  considered  periodic  and  the 
horizontal  boundary  conditions  are  therefore: 


0 u,v,w,p',p')(x+Lx,y,z,t )  , 

(8) 

(' u,v,w,p,p)(x,y,z,t )  = 

(u,v,w,p',p')(x,y  +  Ly,z,t)  ■ 

(9) 

The  bottom  boundary  is  a  solid  wall  with  no  slip  boundary  conditions: 


u(x,y,0,t)  =  0,  v(x,  y,  0,  t)  =  0,  w(x,y,Q,t)=Q 
The  top  boundary  is  a  free-slip  free-surface: 


(10) 


du 

dz 


(■ x,y,Lz,t ) 


=  0  ,  ^ 


dv 

dz 


{ x,y,Lz,t ) 


=  0 


w(x,y,L-,t)  =  0 


(11) 


The  density  perturbation  is  subject  to  a  Dirichlet  boundary  condition  at  both 
vertical  boundaries: 


p(x,  y,  0,  t)  =  p'(x,  y,  Lz,  t)  =  0  .  (12) 

The  boundary  conditions  for  the  pressure  are  of  purely  numerical  nature  and 
their  discussion  is  thus  deferred  to  §3.1. 


3  Numerical  Method 


3.1  Temporal  Discretization 


For  the  temporal  discretization  of  eqs.  (l)-(3),  a  high-accuracy  pressure  projec¬ 
tion  scheme  [38]  is  employed,  which  has  been  well  tested  in  the  recent  spectral 
method  literature  (e.g.  [37,1]).  According  to  this  scheme,  if  one  integrates  eqs. 
(l)-(3)  in  time  from  level  tn  to  tn+\  one  obtains  the  following  semi-discrete 
equations  decomposed  into  three  fractional  steps  for  u: 
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^  Ji  —  1  ' 

u  -  E,_ o  “iu 


(13) 


n—q 

At 


Je~  1 


q= 0 


u  —  u 

At 


-V</>n+1 


(14) 


To11 


n.+  l 


U 


At 


=  z/L(u 


n+ 1  \ 


(15) 


The  splitting  procedure  for  p'  consists  of  two  steps  analogous  to  eqs.  (13)  and 
(15).  In  eqs.  (13)-(15),  a  3rd  order  backward  differentiation  scheme  (BDF3 
with  Ji  =  3)  [9]  is  used  to  discretize  the  temporal  derivative  and  enhance 
the  temporal  accuracy  of  the  splitting  approach.  The  viscous  operator  L  is 
treated  fully  implicitly.  The  weakly  dissipative  nature  of  such  an  approxima¬ 
tion  is  hclfpul  in  stability-sensitive  under-resolved  problems.  The  non-linear 
terms  N  are  advanced  in  time  via  a  third  order  stiffly  stable  scheme  (SS3  with 
Je  =  3)  [38]  allowing  for  maximum  value  of  a  stable  timestep.  The  values  of 
the  coefficients  aq,  (3q  and  70  for  a  BDF3-SS3  scheme  may  be  found  in  Karni- 
adakis  et  al.  [38]. 

The  quantity  (pn+1 : 


£71+1 

I  Vp'dt  =  AtV(pn+1  .  (16) 

in 

is  an  intermediate  scalar  field  that  ensures  that  the  final  velocity  un+1  is 
incompressible.  Note  that  in  eq.  (14),  it  is  assumed  V-u  =  0  and  the  Poisson 
equation  is  solved  for  the  pressure: 


vy+i  =  v.(^)  ,  (17) 

The  boundary  conditions  for  u  (10)-(11)  are  enforced  in  eq.  (15)  and  an  anal¬ 
ogous  approach  is  followed  for  p' .  However,  special  care  must  be  taken  for 
the  correct  choice  of  intermediate  boundary  conditions  for  (j)  in  eq.  (17).  An 
intermediate  boundary  condition  that  ensures  high  temporal  accuracy  for  the 
given  scheme  is  [38]: 
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d(j)n+1 

dz 


b 


(18) 


Je-  1 


n—q i 


9=0 


Je-1 


^  /?Jz/V  x  (V  x  w) 


n—q 


q= o 


where  |&  denotes  z  =  0,  Lz  and  the  coefficients  /3q  have  the  same  value  as  in  the 
SS3  scheme  of  eqs.  (13)-(15).  Guermond  and  Shen  [24]  prove  that  the  above 
splitting  scheme  (eqs.  (13)-(15)  )  and  (18))  is  equivalent  to  the  rotational  form 
of  a  velocity-correction  projection  scheme  whose  second  order  variant  exhibits 
0(At2)  accuracy  in  both  u  and  <t>. 


3.2  Spatial  Discretization 


In  the  periodic  horizontal  direction,  Fourier  spectral  discretization  is  used  with 
Nx  and  Ny  Fourier  modes  in  the  longitudinal  and  spanwise  direction,  respec¬ 
tively.  Horizontal  derivatives  are  calculated  in  a  straightforward  fashion  in 
Fourier  spectral  space.  In  the  vertical  direction,  the  computational  domain  is 
partitioned  into  M  subdomains  of  variable  height  H k  and  order  of  polynomial 
approximation  Nk  {k  —  1,  ...  M)  (figure  2).  Within  each  subdomain,  Legendre 
spectral  discretization  is  used  and  for  the  specific  problem  under  considera¬ 
tion  Nk  is  fixed  and  equal  to  a  fixed  value  N  in  all  subdomains  1 .  In  each 
subdomain,  any  function  f(z)  may  be  approximated  on  the  Gauss  Legendre 
Lobatto  grid  as  [6]: 


Nk- 1 

m  =  £  jjPfa)  •  (w) 

3=0 

where  fj  is  the  amplitude  of  the  j-th  Legendre  mode  and  Pj{z)  the  j-th  order 
Legendre  approximating  polynomial  [6].  Spectral  multidomain  methods  are 
designed  as  collocation  methods  [6]  and  in  this  study  f(z)  is  approximated  in 
nodal,  not  modal,  form  [4,22]: 


Nk- 1 

f(z)  =  f(zj)lj(z)  ■  (2°) 

3=0 

where  lj(z)  is  the  j-th  order  Lagrange  interpolating  polynomial  (cardinal  func¬ 
tion)  [6,4],  The  discrete  vertical  derivative  of  f(z)  is  computed  by  linearly 
mapping  the  global  vertical  coordinate  z  within  a  specific  subdomain  onto  a 

1  In  this  manuscript,  the  symbol  N  may  represent  two  different  quantities,  the 
order  of  polynomial  approximation  in  a  subdomain  or  the  Brunt- Vaisala  frequency. 
In  a  given  section  of  the  manuscript,  the  quantity  N  specifically  refers  to  is  obvious 
from  the  context  of  the  accompanying  text. 
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local  coordinate  (  G  [—1, 1]  with  £  =  —  1  and  1  corresponding  to  the  bottom 
and  top  end-point  of  the  subdomain,  respectively: 


d[_ 

dz 


(*0 


d  d( 


2  Nk~1 

jf  f(zj)Dv 

k  3=0 


(21) 


where  DtJ  is  the  Legendre  spectral  differentiation  matrix.  Conventionally,  D%J 
is  computed  as  a  function  of  Pj  [6].  However,  such  a  calculation  is  error-prone 
especially  for  higher  order  derivatives,  under-resolved  functions  and  high  val¬ 
ues  of  Nk  [14].  Such  errors  may  prove  to  be  catastrophic  in  the  family  of 
splitting  schemes  that  do  not  require  boundary  conditions  for  the  pressure 
[11]  due  to  the  explicit  evaluation  of  third-order  vertical  derivatives.  In  the 
splitting  scheme  described  in  §3.1,  as  a  safeguard  against  errors  induced  by 
spectral  differentiation  an  alternative  vertical  derivative  calculation  technique, 
described  by  Costa  and  Don  [8],  is  used.  The  above  technique  and  the  classi¬ 
cal  pseudospectral  approach  [6]  are  used  to  compute  all  the  derivatives  in  the 
estimation  of  the  non-linear  term  of  eq.  (13). 


In  any  case,  as  discussed  in  §1.3,  at  high  Re,  the  existing  splitting  scheme  and 
spectral  multidomain  discretization  is  prone  to  significant  numerical  instabil¬ 
ities  most  evident  at  the  boundaries  and  subdomain  interfaces.  The  source  of 
these  instabilities  are  aliasing  effects  driven  by  thin  physical/numerical  bound¬ 
ary  layers.  The  instabilities  generate  spurious  energy  with  increasingly  higher 
and  higher  frequency  content  which  has  a  catastrophic  effect  on  the  long-term 
integration  of  the  governing  equations  [22],  It  is  imperative  to  apply  stabiliza¬ 
tion  methodologies,  namely  penalty  methods  and  in  addition,  spectral  filtering 
and  strong  adaptive  interfacial  averaging. 


3.3  Spectral  Multidomain  Penalty  Methods 


The  recent  thrust  in  development  of  penalty  methods  originated  from  the 
need  to  implement  appropriate  boundary  conditions  for  single  domain  wave- 
dominated  dissipative  problems  in  a  manner  that  is  easy,  stable  and  accurate 
[27].  The  direct  enforcement  of  the  boundary  conditions  does  not  guarantee 
that  the  equation  is  satisfied  arbitrarily  close  to  the  boundary.  At  high  Re, 
thin  boundary  layers  and  consequently,  unresolved  gradients  develop  which 
contaminate  the  flow  with  Gibbs  oscillations  through  aliasing.  As  already 
discussed  in  §1.4,  spectral  filtering  does  damp  these  oscillations  in  the  do¬ 
main/subdomain  interior  but  is  least  effective  near  the  boundaries/interfaces. 
When  using  a  collocation  method  as  is  often  done  with  Legendre/Chebyshev 
schemes,  this  problematic  behavior  may  be  overcome  by  collocating  both  the 
equation  and  the  boundary  condition  at  the  physical  boundary,  the  latter  mul- 
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tiplied  by  a  penalty  coefficient.  As  a  result,  the  eigenvalues  of  the  associated 
differential  operators  are  modified  [27]  and  a  smooth  numerical  solution  is  ob¬ 
tained  subject  to  an  error  (penalty)  equal  to  the  order  of  the  spectral  scheme. 
The  spectral  accuracy  of  the  numerical  scheme  is  negligibly  impacted  and  one 
may  compute  stably  the  high  Re  internal  dynamics  of  the  flow  without  having 
to  resolve  the  thin  physical/numerical  boundary  layers.  Consequently,  penalty 
methods  are  a  key  enabling  technique  for  the  simulation  of  localized  high  Re 
turbulence. 

In  the  case  of  spectral  multidomain  discretization,  the  core  of  the  penalty 
method  is  similar  to  that  of  the  single  domain  case  [28].  Each  individual  sub- 
domain  may  be  treated  as  a  single  domain  whose  interfaces  support  patching 
conditions  that  act  as  open  boundary  conditions  through  which  information  is 
exchanged  with  adjacent  subdomains.  The  exact  nature  of  the  patching  con¬ 
ditions  depends  on  whether  the  equation  is  parabolic  or  hyperbolic.  In  both 
cases,  the  patching  conditions  are  such  that  C0  and  C\  continuity  is  enforced 
only  weakly  since  the  interface  of  two  adjacent  subdomains  corresponds  to 
the  same  location  in  physical  space  but  to  two  separate  grid-points.  Thus, 
the  penalty  method  is  inherently  discontinuous  and  it  is  this  weak  continuity 
on  the  numerical  grid  that  allows  for  a  more  stable  and  smooth  numerical 
solution  at  the  subdomain  interfaces. 

In  terms  of  the  splitting  scheme  outlined  in  §3.1,  the  penalty  method  is  ap¬ 
plied  at  two  different  levels  in  the  incompressible  Navier-Stokes  equations. 
The  explicit  non-linear  term  advancement  is  treated  as  a  hyperbolic  equation 
whereas  the  implicit  viscous  term  treatment  as  a  parabolic  equation  (in  this 
section  all  subsequent  equations  are  written  as  a  function  of  the  u-velocity 
without  loss  of  generality): 

^  =  N(«)  ,  (22) 


Qi  =  v L(«)  .  (23) 

The  temporal  derivatives  in  the  equations  (22)- (23)  are  only  approximations 
to  those  appearing  in  eqs.  (13)-(15). 

Consider  first  eq.  (22).  Here  Gibbs  oscillations  at  subdomain  interfaces  must 
be  suppressed.  The  boundary  conditions,  imposed  strictly  in  eq.  (15),  are 
not  incorporated  in  the  penalty  treatment  of  the  hyperbolic  problem.  Well- 
posedness  dictates  that  in  each  subdomain  of  index-/;;  and  uniform  order 
Nk  =  N  (all  subsequent  penalty  formulations  are  valid  for  variable  as 
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well),  eq.  (22)  be  subject  to  the  following  set  of  patching  conditions  at  the 
bottom  and  top  subdomain  interfaces,  respectively: 


auQ=g^(t)  ,  where  gk(t)  =  aukNl,  (24) 

7 ukN  =  gk2(t)  ,  where  g%(t)  =  yu k0+1. 

The  penalty  formulation  of  (22),  in  physical  space,  is: 


duk 

~dt 


N  {uk) 


t?Q 


fc/n- 1 

k  ' 


au0 


T2  Qk(Z; 


7  ukN 


,  (25) 


where 


Qk(zi)  =  ho,  Qt(zk)  =  SlN  .  (26) 

here  Sij  is  the  Kronecker  delta  function  with  subscript  i  corresponding  to 
the  collocation  point  zk.  Following  the  treatment  by  Hesthaven  of  patching 
conditions  as  localized  open  boundary  conditions  (locally,  each  interface  expe¬ 
riences  inflow  or  outflow),  the  values  of  the  coefficients  a,  7,  rk  and  rk  and  are 
determined  by  considering  the  value  of  the  vertical  interfacial  velocities  Wq 
and  W%  at  the  previous  timesteps.  In  the  discretized  form  of  the  equations 
(where  Q  -see  §3.2-  is  the  vertical  coordinate),  if  n  is  the  vector  normal  to  the 
subdomain  interface,  the  following  cases  are  distinguished: 


Wq  ■  h0  >  0  : 

a 

=  0, 

Tk 

=  0  , 

Wq  •  h0  <  0  : 

a 

=  Wo,. 

Tk 

—  1  2 

2t 0  Hk 

W^j  ■  hjv  >  0  : 

'■  7 

=  0, 

t2 

=  0  , 

Wx  ■ hN  < 0  : 

:  7 

=  \wN\ 

_  1  2 
2uj  Hk 

where: 


2 

N(N  +  1) 


(27) 


and  Hi ,  is  the  height  of  the  k-th  subdomain.  The  values  of  rk  and  Tk  are  what 
is  used  in  the  discretized  formulation  of  (25).  Eqs.  (24)  and  (25)  imply  weak 
Co  continuity  at  the  subdomain  interfaces  due  to  the  normal  advective  flux  of 
momentum  /  mass. 


The  penalty  formulation  of  (23)  is  significantly  different  from  that  of  (22)  due 
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to  its  parabolic  nature.  Eq.  (23)  is  solved  in  Fourier  space  and  its  discretized 
form  for  each  pair  of  horizontal  wavenumbers  (kx,  ky)  is 


D2u  -  ( k2x  +  k2y 


7o  x  u 

—— )u  = - — 

vAt 7  vA  t 


where 


(28) 


and  u  is  the  value  of  the  velocity  at  timestep  n  +  1 .  It  is  immediately  obvious 
that  in  the  case  of  significantly  high  Re,  (28)  has  excessively  thin  numerical 
boundary  layers  of  thickness  0( \/ vA t)  whose  resolution  requires  prohibitively 
high  number  of  grid  points  at  the  boundary  to  be  resolved  [4,9].  As  a  result,  a 
singularly  perturbed  boundary  value  problem  arises  which  has  serious  reper¬ 
cussions  on  the  stability  of  the  numerical  solution.  In  addition,  at  high  Re, 
strict  enforcement  of  Co  and  C\  continuity  in  the  solution  of  (28)  over-excites 
the  higher  Legendre  modes  resulting  in  spurious  oscillatory  behavior  at  the 
subdomain  interfaces.  For  the  case  of  the  under-resolved  Navier-Stokes  equa¬ 
tions,  the  neglect  of  any  penalty  terms  for  the  non-linear  term  computation, 
compounds  these  oscillations  through  aliasing  effects. 

To  derive  the  penalty  form  of  (28)  one  defines  a  small  quantity: 


e  =  (K  +  ky 


7o  \  — i 
vAt 


and  recasts  (28)  as: 


(29) 


u 

u  —  eD2u  =  e — —  =  F  .  (30) 

uAt  K  ’ 

The  boundary  operators  in  each  subdomain  are  determined  by  well  posedness 
and  the  elliptical  nature  of  (30)  and  are  given  by 


au0  -  =  01 W  ’  where  9i(t)  =  aMv1  -  ^e~gZk-[  >  (31) 

JutN  +  SeiA-  =  SUt)  .  "here  9*(i)  =  7'C‘  +  ■ 

Thus,  the  penalty  formulation  of  (30)  is 
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eD2u  —  u  +  F 


(32) 


FkQ~k  I 


z ■ 


aun 


-  T. 


:Qt 


r)i  i  ^ 

Be ^ 
P  dzk 

du\ 


N 


iuk  +  5e  dzk 


9i(t) 


=  0 


where  Q^(zk)  and  Q^(zk)  are  defined  in  (26).  It  is  obvious  that  (31)  and 
(32)  enforce  weak  Co  and  C'\  continuity  at  the  subdomain  interfaces.  At  the 
physical  boundaries,  g\{t)  =  gi(t)  and  =  g2 (t)  where  gi(t)  and  g2(t)  are 
the  boundary  conditions  (10)-(12).  In  the  case  of  non- homogeneous  boundary 
conditions  (as  may  happen  with  the  pressure,  discussed  below),  gi(t)  and 
g2(t)  are  not  the  exact  physical  conditions  but  a  modified  form  that  ensures 
the  correct  term  balance  in  (31).  The  coefficients  a,  /3,  7  and  5  are  set  to  1 
at  each  subdomain  interface  and  at  the  physical  boundaries,  their  values  are 
set  to  satisfy  the  physical  boundary  conditions.  The  discretized  form  of  (32) 
supplemented  by  (31)  results  in  the  solution  of  a  linear  system  of  equations 
where  the  coefficient  matrix  is  band-diagonal  with  bandwidth  N  and  may  be 
solved  through  a  variety  of  existing  fast  LU  solvers.  Finally,  the  coefficients  rf 
and  r|  for  the  discrete  equations  (where  (  is  the  mapped  vertical  coordinate) 
are  determined  by  specific  stability  constraints,  for  the  following  cases  [27]: 


(1)  Physical  Boundary  with  Dirichlet  Condition: 


Ti  = 


au) 


2 

2^H, 


r,M  = 


~(- 
5lo2  V  H, 


M 


where  u  is  defined  in  (27). 

(2)  Physical  Boundary  with  Neumann  Condition: 


1  2 


Ti  = 


rM  - 
'2 


I  2 


Ho  z  5u  Hm 

(3)  Subdomain  Interface  with  Robin  Patching  Condition: 


ri  = 


tue/3  1 

1 


To  = 


cueS  L 


6  H-  2^1  —  2  'S^J  H-  €K,\ 
e  +  2k2  —  2  \ /k  2  +  ek2 


2 

j H~k 
2 


(33) 


(34) 


(35) 


where  /sq  =  Lva//3  and  n2  =  LU'f/S .  In  the  numerical  implementation  of  a 
Neumann  boundary  condition,  the  value  indicated  by  (34)  should  be  strictly 
adhered  to.  For  Dirichlet/  Robin  conditions,  the  values  indicated  in  (33)  and 
(35)  are  the  lower  bound  of  associated  stability  requirements.  In  the  simula¬ 
tions  of  this  paper,  the  penalty  coefficients  for  a  Dirichlet  condition  are  kept 
fixed  at  the  value  prescribed  above.  However,  the  upper  stability  bound  on 
the  penalty  coefficient  of  a  Robin  patching  condition  is  of  0(e-1)  or  greater 
than  the  corresponding  lower  bound  [27].  It  was  found  that  using  a  value 
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equal  to  that  of  (35)  multiplied  by  a  factor  of  O^e-1)  increased  the  stability 
of  the  scheme  at  the  subdomain  interfaces.  An  increased  value  of  an  inter¬ 
facial  penalty  coefficient  corresponds  to  stronger  enforcement  of  the  patching 
conditions  and  thereby  more  local  dissipation  leading  to  increased  stability. 

The  penalty  treatment  of  (28)  is  applied  to  the  Poisson  equation  for  the  pres¬ 
sure  (17)  in  spectral  space  for  any  given  pair  of  horizontal  Fourier  wavenum¬ 
bers  ( kx ,  ky)  because  (17)  is  a  one-dimensional  Helmholtz  equation.  The  only 
exception  is  the  (0,  0)  Fourier  mode  of  0  which,  however,  is  irrelevant  to  the 
computation  because  the  corresponding  mode  of  w  is  zero  due  to  incompress¬ 
ibility  and  the  boundary  conditions  for  w  while  the  (0,  0)  modes  of  u  and  v 
are  not  0-dependent. 

Note  that  the  use  of  SMPM  discretization  allows  maximum  flexibility  in  the 
choice  of  boundary  operators  (Dirichlet/ Neumann/ Robin  and  homogeneous/non- 
homogeneous) .  Thus,  the  existing  spectral  multidomain  model  is  amenable  to 
the  study  of  geophysical  flows  with  complex  boundary  forcings  (wind  stress/buoyancy 
flux) . 


3-4  Spectral  Filtering 


For  the  grid  resolutions  considered  in  this  paper,  penalty  methods  allow  an 
increase  of  the  value  of  Re  by  roughly  two  orders  of  magnitude.  However,  even 
such  an  increase  is  not  sufficient  to  attain  values  of  Re  typical  of  laboratory 
wake  experiments  and  beyond.  Some  form  of  a  hyperviscous  operator  is  needed 
to  ensure  stability  for  these  high  Re.  In  practice  what  is  used  instead,  which 
has  the  same  effect  but  avoids  additional  stiffness  and  subsequent  timestep 
limitations,  is  the  application  of  a  low-pass  spectral  filter  on  the  numerical 
solution  whose  convergence  rate  is  thereby  enhanced  [22],  As  long  as  certain 
basic  requirements  [41,22]  are  satisfied,  no  unique  choice  of  a  filter  function  is 
required  and  in  this  study,  an  exponential  filter  [22]  is  used: 


|  1,  0  <  k  <  kc  , 


a(k)  =  { 


exp 


—  a 


(  k-kc  \ 

(.iV-fccJ 


p  i 


kc  <  k  <  N 


(36) 


where  p  is  the  filter  order,  kc  the  filter  lag  and  a  =  — Iu£m  with  eM  being  the 
machine  precision.  The  filtered  solution  f*  may  now  be  expressed  in  terms  of 
the  modes  of  the  numerical  solution  as: 
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(37) 


Nk- 1 

=  5Z  a(ki)fiPj(zi) 

3=0 

where  kj  is  the  j-th  discrete  Legendre  mode.  An  analogous  expression  may  be 
written  for  filtering  in  Fourier  space  [22], 

A  common  concern  with  the  implementation  of  spectral  filtering  in  spatially 
continuous  spectral  element  methodologies  is  that  filtering  does  not  preserve 
the  patching  and  boundary  conditions  and  thus  specific  measures  need  to  be 
adopted  [41,3,1].  Such  a  concern  does  not  exist  when  using  the  inherently 
discontinuous  penalty  method  because  the  error  induced  by  the  filtering  op¬ 
eration  is  of  the  same  order  as  the  penalty  scheme  [28],  i.e.  minimal. 

In  the  incompressible  spectral  multidomain  solver  presented  in  this  paper, 
spectral  filtering  is  applied  at  three  levels  when  advancing  the  solution  from 
time  level  (n)  to  level  (n  +  1).  First,  to  eliminate  aliasing  effects,  filtering  is 
applied  in  both  Fourier  and  Legendre  space  after  advancing  the  non-linear 
terms  in  (13).  The  solution  to  the  Poisson  equation  (17)  is  filtered  in  the  ver¬ 
tical  direction  to  smooth  out  any  errors  induced  at  the  subdomain  interfaces 
due  to  the  discrete  estimate  of  dw/dz.  Finally,  the  solution  of  (15)  is  filtered 
in  the  vertical.  The  order  of  the  Legendre  filter  is  the  same  in  all  three  levels. 

For  the  Re  under  consideration,  use  of  spectral  filtering  is  the  inevitable  price 
when  confronted  with  the  associated  high-degree  of  under-resolution.  Filtering 
should  not  be  viewed  as  a  waste  of  resolution  because  a  significant  percent¬ 
age  of  modes  of  the  numerical  solution  (0(40  —  50%)  in  Legendre  space  and 
0(70  —  80%)  in  Fourier  space)  are  unaffected  by  the  filtering  procedure  thus 
preserving  the  high  order  of  accuracy  of  the  spectral  scheme.  The  filtered  up¬ 
per  band  of  modes  is  necessary  to  ensure  stability  and  effect  energy  transfer 
to  the  unresolved  scales. 


3.5  Strong  Adaptive  Interfacial  Averaging 


Although  filtering  enhances  even  further  the  stability  properties  of  a  spectral 
multidomain  penalty  scheme,  the  subdomain  interfaces  occasionally  exhibit  a 
propensity  for  numerical  instability.  Within  each  subdomain,  using  the  trans¬ 
formed  coordinate  application  of  the  exponential  filter  of  (36)  is  equivalent 
to  incorporation  in  the  equation  of  the  differential  operator  [22]: 


dw5 


d 


acJ 


(38) 
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Obviously,  (38)  indicates  that  spectral  filtering  is  least  effective  near  the  sub- 
domain  interfaces  where  the  highest  Legendre  modes  are  inherently  most  os¬ 
cillatory  [6].  To  prevent  the  development  and  subsequent  growth  of  any  such 
interfacial  singularities  a  simple  adaptive  interfacial  averaging  technique  is 
used.  Once  the  final  value  of  the  numerical  solution  for  (u,  v,  w,  p')  has  been 
obtained  in  physical  space,  at  the  top  interface  of  the  k-th  subdomain,  when 
the  following  criterion  is  met  [13]: 


u, 


fc+i 


—  u 


AM 


U  i 


k+ 1 


+  U 


>Ca 


N I 


an  averaging  operation  is  performed: 


(39) 


<4  =  4+1  =  )(d+1 +>*_,)  .  (40) 

The  coefficient  Cave  is  set  equal  to  0.005.  Equation  (40)  was  found  to  be  more 
effective  than  the  technique  used  by  Don  et  al.  [13].  Including  the  interfacial 
points  in  the  averaging  and  even  more  so  restricting  the  averaging  to  the 
interface  was  found  to  not  eliminate  troublesome  spikes  in  the  solution.  The 
averaging  essentially  adds  some  weak  local  artificial  dissipation  to  the  scheme 
and  is  usually  required  for  only  a  small  fraction  (0(0.1%))  of  the  subdomain 
interfaces. 


4  Results 


4-1  Problem  Geometry 


The  base  flow  chosen  to  validate  the  spectral  multidomain  penalty  model  is 
a  momentumless  stratified  turbulent  wake.  Such  a  flow  corresponds  to  the 
mid-to-late  time  wake  of  a  sphere  of  diameter  D  towed  with  a  velocity  U  in  a 
linear  density  stratification  of  frequency  N ,  where 


TV2  = 

Po  dz 

Note  that  the  spatial  discretization  does  not  account  for  the  sphere  and  fo¬ 
cuses  only  on  the  flow  generated  in  its  wake.  The  computational  domain  of 
dimensions  Lx  x  Ly  x  Lz  is  shown  in  figure  1  and  corresponds  to  a  window 
fixed  in  space  with  respect  to  the  moving  sphere,  as  do  the  laboratory  mea¬ 
surements  [50].  Behind  the  sphere,  the  wake  is  considered  to  be  statistically 
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stationary.  The  periodicity  assumption  in  the  x-direction  is  valid  because  the 
length  of  the  computational  domain  is  much  smaller  than  the  total  wake  length 
[12].  The  spanwise  periodicity  assumption  is  also  valid  provided  the  horizontal 
lengthscale  of  the  wake  does  not  become  excessively  large.  Due  to  spanwise 
periodicity,  internal  waves  radiated  by  the  wake  re-enter  the  computational  do¬ 
main  and  thus  after  a  certain  point  later  in  time,  the  internal  wave  held  does 
not  correspond  to  its  experimental  counterpart.  However,  by  that  time  the 
how  is  efficiently  decomposed  into  internal  waves  and  quasi-two-dimensional 
vortical  modes  and  analysis  of  the  latter  [44]  is  possible  in  isolation. 


4-2  Initialization:  Replacing  the  Sphere 


Although  the  sphere  is  not  accounted  for  in  the  computation,  its  effect  must  be 
incorporated  in  the  initial  condition.  The  initial  how  held  is  the  superposition 
of  an  axisymmetric  Gaussian  mean  velocity  prohle  and  a  turbulent  fluctuation 
held: 


u (x,y,z,t)  =  Ux(y,z,t)  +  u'(x,y,z,t)  .  (42) 

The  X  subscript  indicates  averaging  in  the  streamwise  direction.  The  mag¬ 
nitude  and  distribution  of  mean  and  fluctuating  velocity  fields  are  specihed 
based  on  available  vertical  prohle  measurements  of  stratihed  sphere  wakes 
[50,47-49]  at  Nt  =  3  (downstream  distance  from  the  sphere,  x/D  =  6). 

The  mean  velocity  prohle  is  a  Gaussian: 


ux(y,z,t) 


U0  exp 


1  (y~\Lv\2 

2  V  LY  ) 


1 (Z-\Lz\2- 
2  V  Lz  )  . 


(43) 


Uq  is  the  maximum  centerline  velocity  and  Ly  and  Lz  the  initial  horizontal 
and  vertical  lengthscales,  respectively.  Initially,  one  sets  Vx  =  Wx  =  0. 


The  fluctuating  prohle  is  also  taken  from  laboratory  data  at  Re  =  5  x  103 
and  Fr  =  4.  The  x-averaged  r.m.s.  (root  mean  square)  distribution  of  the 
fluctuating  velocity  is  assumed  to  be  axisymmetric  and  equipartitioned  among 
its  three  components: 


%(h  t)  =  v'x  =  w'x 


U0 


(44) 


where  r  =  {{y  —  y0 )2  +  (z  —  ^o)2)1^2-  Although  the  assumptions  of  axisymme- 
try  and  energy  equipartition  are  debatable  at  Nt  =  3,  horizontal  prohles  from 
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experiment  are  available  from  Nt  =  9  onwards  and  a  small  degree  of  extrapo¬ 
lation  is  required  in  specifying  the  initial  condition  in  the  horizontal  direction 
in  the  simulations  (see  §4.4).  A  delayed  onset  of  asymmetry  may  result  but 
its  effect  on  later  wake  statistics  is  negligible.  Ideally,  one  would  initialize  the 
model  with  runs  from  simulations  of  stratified  flow  around  a  sphere  but  a 
spectral  multidomain  scheme  for  such  a  simulation  would  be  a  challenge  of  its 
own.  In  addition,  low-order  scheme-  based  studies  of  sphere  wakes  have  been 
restricted  to  Re  =  200  [26]. 

The  three-dimensional  fluctuating  velocity  held  is  constructed  as  spectrally 
random  noise  in  three-dimensional  Fourier  space  with  a  k~5^3  energy  spec¬ 
trum.  An  inverse  Fourier  transform  is  applied  to  convert  the  noise  into  physical 
space  and  in  the  vertical  the  fields  are  projected  on  the  non-uniform  Gauss- 
Lobatto-Legendre  grid  of  each  subdomain.  Finally,  the  data  is  windowed  in 
the  envelope  of  the  r.m.s.  profile  of  eq.  (44).  Use  of  white  noise  is  avoided 
because  it  is  unphysical  and  detrimental  to  the  stability  of  the  numerical  so¬ 
lution. 

The  initial  density  held  is  assumed  to  be  an  unperturbed  linear  density  strati¬ 
fication.  Transients  due  to  the  gravity-induced  adjustment  of  the  density  held 
to  the  turbulence  [10]  have  no  significant  effect  on  later  wake  dynamics. 

Initially,  the  fluctuating  and  mean  velocity  helds  are  uncorrelated.  Thus,  a 
preliminary  “relaxation”  simulation  [12]  is  run  to  generate  a  physically  re¬ 
alistic  velocity  held.  During  relaxation,  the  how  is  forced  to  maintain  con¬ 
stant  mean  and  r.m.s.  huctuating  velocity  prohles  according  to  eqs.  (43)  and 
(44),  while  the  spatial  distribution  of  the  turbulent  fluctuations,  and  thus  the 
Reynolds  stresses,  is  allowed  to  vary. 


4- 3  Run  Description 


The  spectral  multidomain  penalty  method  model  was  used  to  simulate  mid-to- 
late  time  momentumless  stratified  turbulent  wakes  at  Reynolds  number  Re  = 
UD/u  —  5  x  103  and  2  x  104  and  internal  Froude  number,  Fr  =  U /( DN )  =  4 
and  oo.  U  and  D  are  the  tow  velocity  and  diameter  of  the  virtual  sphere  to 
which  the  initial  wake  condition  at  Nt  =  3  would  correspond  to.  Re  is  mod¬ 
ified  by  changing  the  value  of  the  molecular  viscosity  v.  Re  =  5  x  103  and 
Fr  =  4  are  typical  values  of  laboratory  experiments  performed  by  Speckling 
and  coworkers  [50,47-49].  The  choice  of  Re  and  Fr  values  for  the  simulations 
was  made  to  illustrate  the  effect  of  the  turbulent  scale  separation  (increasing 
with  Re)  and  buoyancy  on  the  dynamics  of  a  turbulent  wake.  Re  =  2  x  104  is 
a  maximum  Re  for  the  resolution  employed  in  this  study,  above  which  under¬ 
resolution  effects  affect  the  robustness  of  the  simulation. 
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The  computational  domain  has  a  horizontal  dimension  of  LxxLy  =  16  D  x  16  D 
and  corresponds  to  a  virtual  stratified  water  tank  of  height  H  =  12 D.  Such 
domain  dimensions  ensure  that  the  domain  exhibits  adequate  length,  width 
and  height  to  allow  for  multiple  streamwise  wavelengths  of  a  vortex  shedding 
instability  and  also  minimize  effects  of  confinement  and  interaction  with  the 
wake’s  periodic  image.  A  uniform  spatial  grid  is  used  in  the  horizontal  direc¬ 
tion  whereas  in  the  vertical  the  spectral  multidomain  discretization  of  figure 
2  is  employed  with  M  =  5  non-uniform  height  subdomains  of  fixed  order 
of  approximation  Nk  =  N.  Increased  resolution  is  available  at  the  energetic 
core  of  the  wake  whereas  less  is  utilized  in  the  less  active,  internal- wave  domi¬ 
nated  ambient.  The  non-uniform  Gauss  Lobatto  Legendre  grid  of  the  topmost 
subdomain  allows  for  enhanced  resolution  of  the  subsurface  region.  The  ra¬ 
tio  of  adjacent  subdomain  heights  Hk/Hk- 1  is  maintained  within  the  interval 
[0.5,2],  Otherwise,  reflections  are  observed  at  the  interfaces.  The  resolution 
used  for  the  Re  =  5  x  103  and  Re  =  2  x  104  runs  is  128  x  128  x  165  and 
128  x  128  x  245,  respectively.  The  low  Re  Fourier  and  Legendre  spectral  fil¬ 
ters  are  ( Pf,Pl )  =  (20,8)  and  at  the  high  Re,  ( Pf,Pl )  =  (16,6)  is  used.  In 
the  higher  Re  run,  increased  under-resolution  is  countered  by  increasing  N , 
while  keeping  Hk  fixed  (p-refinement),  and  reducing  filter  order.  The  adequacy 
of  resolution  for  the  Re  =  5  x  103  and  2  x  104  is  evident  in  figure  3  which 
shows  one-dimensional  spectra  of  the  turbulent  kinetic  energy  at  three  differ¬ 
ent  times  for  the  Fr  =  oo  simulations.  The  energy  spectra  drop  off  smoothly 
at  the  smallest  resolved  scales,  where  no  sign  of  artificial  energy  accumula¬ 
tion,  an  unwanted  feature  of  under- resolution,  is  observed  [4],  Despite  very 
similar  initial  spectra  at  t  =  Os  (a  possible  effect  of  the  relaxation  procedure), 
Re  =  2  x  104  develops  different  features,  primarily  the  presence  of  a  more 
extended  inertial  range  of  nearly  a  decade  in  span. 

The  computational  timestep  At  was  chosen  as  such  that  the  CFL  stability 
criterion  be  obeyed  in  all  three  spatial  directions  for  a  3rd  order  stiffly  stable 
scheme.  The  following  requirements  are  imposed: 


^Urnax  <  ^Vrnax  <  ^ 

Ax  Ay 

In  all  runs,  at  t  —  0,  A t/(D/U)  =  0.04  is  used.  At  later  times  of  a  simulation, 
when  the  CFL  constraint  is  relaxed,  the  run  may  be  restarted  with  a  higher 
timestep.  It  is  imperative  that  3rd  (not  lower)  order  BDF  be  used  in  the  restart 
scheme.  In  the  temporally  under-resolved  simulations  of  interest,  lower  order 
BDF  approximation  of  the  temporal  derivative  produces  excessive  artificial 
dissipation  and  non-physical  decrease  of  the  kinetic  energy  of  the  flow.  During 
the  first  two  timesteps  of  the  restart,  variable  timestep  BDF3  [25]  and  3rd 
order  Adams-Bashforth  schemes  [46]  are  used. 


w 

~Az 


<  0.6 


(45) 
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4-4  Setting  the  Initial  Wake  Lengthscales  and  Velocities 


The  initial  values  of  the  coefficients  in  the  mean  and  fluctuating  velocity  pro¬ 
file  eqs.  (43)  and  (44)  are  either  taken  directly  or  extrapolated  from  the  data 
of  Spedding  [47-49] .  Given  Spedding’s  observation  that  the  horizontal  growth 
rate  of  a  stratified  wake  is  equal  to  that  of  its  unstratified  counterpart,  one 
may  use  later  time  data  for  Fr  =  4  and  at  x/D  —  6,  set  Ly / D  =  0.4.  The  ini¬ 
tial  time  Nt  =  3  falls  within  the  non-equilibrium  (NEQ)  regime  during  which 
Uq  follows  a  f-0'25  power  law.  Thus,  one  can  estimate  from  Spedding’s  data 
available  at  Nt  =  5,  a  value  of  Uq/U  =  0.1479  at  Nt  =  3.  Using  laboratory 
data  at  Nt  =  5  and  the  assumption  that  the  constant  value  of  Lz  during  the 
NEQ  regime  extends  as  far  back  as  Nt  =  2,  as  observed  in  the  self-propelled 
body  experiments  of  Lin  and  Pao  [42],  one  estimates  Lz  =  0.4  =  Ly  at 
Nt  =  3.  Axisymmetry  in  mean  wake  geometry  is  plausible  during  the  earlier 
phases  of  the  buoyancy-influenced  NEQ  regime.  In  addition,  it  motivates  the 
choice  of  an  axisymmetric  fluctuating  velocity  profile  in  §4.2.  Laboratory  data 
from  vertical  slices  indicate  that  at  t  —  0  in  (44),  ay  =  0.03375,  /3\  =  1/15, 
rg/D  =  0.35  and  rp/D  =  0.1.  The  Fr  =  oo  simulations  utilize  the  same  initial 
condition.  It  is  questionable  whether  two  towed  body  wakes  at  Fr  =  4  and 
Fr  =  oo  at  identical  Re  and  Pr  have  the  same  characteristic  lengths  and 
velocities  at  x/D  —  6.  However,  the  focus  of  this  paper  is  the  reproduction  of 
the  stratified  wake  experimental  data  and  the  sole  concern  with  the  Fr  =  oo 
wakes  is  that  they  obey  the  appropriate  temporal  power  laws. 

All  simulations  are  run  until  the  dimensions  of  the  characteristic  vortical 
structures  (see  §4.5)  become  large  in  comparison  to  those  of  the  computa¬ 
tional  domain.  The  spanwise  dimension  of  the  domain  is  most  restrictive  in 
the  Fr  =  4  runs,  whereas  the  Fr  =  oo  simulations  experience  confinement 
primarily  clue  to  the  height  of  the  domain. 


4-  5  Model  Validation:  Qualitative  Results 


Figure  4  shows  isosurfaces  of  the  vertical  vorticity  uz  for  the  Fr  =  4  simu¬ 
lations  at  time  Nt  =  56  (x/D  =  112).  The  formation  of  large  horizontal  (“ 
pancake”)  eddies,  three-dimensional  structures  with  a  large  aspect  ratio,  typ¬ 
ical  of  the  late  evolution  of  stratified  turbulence  [16,50,44]  is  immediately  ev¬ 
ident.  Pairings  between  pancakes  of  like-signed  vorticity  gradually  take  place 
as  the  eddies  grow  in  size  but  diminish  in  population.  At  higher  Re,  a  slight 
distortion  in  pancake  geometry  is  observed,  although  the  “pancakes”  exhibit  a 
larger  aspect  ratio  [16].  Fr  =  oo  isosurfaces  of  vorticity  magnitude  (figure  5) 
exhibit  much  less  organized  structure  of  a  smaller  characteristic  lengthseale. 
The  vorticity  field  is  organized  in  randomly  oriented  vortex  tubes.  Thinner 
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and  shorter  tubes  characterize  the  Re  =  2  x  104  unstratified  visualizations  as 
an  increasingly  finer  structure  is  visible  [12].  As  expected  in  an  unstratified 
turbulent  flow,  the  features  of  coherent  structures  in  the  vorticity  field  are 
highly  dependent  on  the  vorticity  magnitude  threshold  [31].  Fr  =  oo  visual¬ 
izations  at  even  lower  cj  thresholds  do  display  quite  a  disorganized  vorticity 
field  but  also  events  comparable  to  vortex  rings  as  shown  in  Gourlay  et  al.  [23] . 

A  y  —  z  section  of  the  horizontal  velocity  divergence, 


A, 


du  dv 
dx  dy 


(46) 


is  shown  in  figure  6  .  At  late  times  and  the  limit  of  low  Froude  number,  Az 
represents  internal  wave  motions  in  the  flow.  At  earlier  times  in  the  develop¬ 
ment  of  the  stratified  wake,  where  separation  of  internal  waves  and  turbulence 
is  not  possible,  A^  is  simply  an  approximate  and  not  exact  descriptor  of  in¬ 
ternal  wave  activity  [50].  The  radiation  of  internal  waves  from  the  stratified 
wake  at  relatively  early  times  (Nt  =  15)  is  evident  in  figure  6.  Although  the 
structure  of  A,  in  the  core  of  the  wake  is  rather  incoherent,  emitted  internal 
wave  rays  are  seen  propagating  away  at  an  angle  of  approximately  45°  from 
the  vertical,  as  in  the  laboratory  data  [49].  Note  that  these  propagating  waves 
experience  no  reflections  at  the  subdomain  interfaces. 


4-6  Model  Validation:  Quantitative  Results 


A  quantitative  validation  of  the  spectral  multidomain  model  is  provided  by 
comparing  results  for  mean  flow  quantities,  the  characteristic  velocity  and 
lengthscales  in  (43),  with  experimental  data  obtained  over  a  range  of  Re  and 
Fr.  The  values  of  Uq,  Ly-  and  Lz  are  obtained  by  performing  a  Gaussian 
fit  according  to  (43)  on  the  x-averaged  mean  profile  calculated  at  individual 
sampling  times  in  the  run  evolution.  The  model  results  are  also  compared  to 
DNS  results  of  Gourlay  et  al.  [23]  at  Re  =  104,  Fr  =  10  and  LES  data  of 
Dommermuth  et  al.  [12]  at  Re  =  104,  Fr  =  4  mixed  model  LES  (their  corre¬ 
sponding  Re  =  105  results  are  not  different).  Comparison  with  results  from 
the  laboratory  and  other  numerical  simulations  is  focused  on  the  stratified 
case.  The  unstratified  runs  are  not  initialized  with  conditions  equivalent  to 
Fr  =  oo  laboratory  wakes  at  x/D  =  6.  Therefore,  validation  of  the  unstrati¬ 
fied  results  is  deemed  sufficient  if  the  timeseries  conform  to  power  laws  derived 
from  similarity  analysis  [51]. 

First,  the  maximum  centerline  velocity  Uq  is  considered  in  figure  7.  Figure  7 
compares  the  spectral  multidomain  results  with  experimental  data  and  other 
DNS/LES  results.  One  of  the  distinguishing  characteristics  of  stratified  wakes 
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is  the  presence  of  low  mean  kinetic  energy  decay  rates  in  the  non-equilibrium 
(NEQ)  regime  for  Nt  €  [2,  50]  [47].  The  SMPM  simulations  produce  the  same 
result.  Both  NEQ  and  quasi  two-dimensional  (Q2D)  power  law  exponents  for 
the  decay  of  Uq  are  within  experimental  uncertainty.  The  DNS  and  mixed 
model  LES  data  exhibit  similar  behavior.  No  significant  change  is  observed 
with  increasing  Re  in  the  spectral  multidomain  data.  The  same  plot  (figure 
7)  also  shows  the  evolution  of  the  SMPM  data  for  Uq  in  the  unstratified  runs. 
For  both  values  of  Re,  the  unstratified  wake  does  decay  as  a  function  of  x~ 2/3 
determined  by  similarity  analysis  [51].  As  also  found  in  the  laboratory  [50,47], 
the  stratified  wake  exhibits  a  higher  value  of  Uq  and  thus  a  higher  mean  energy 
density  due  to  the  constraining  of  the  vertical  wake  scale  Lz  by  the  stratifi¬ 
cation  in  the  NEQ  regime  (see  below). 

The  evolution  of  the  wake  horizontal  lengthscale  Ly  is  shown  in  figure  8a. 
Experimental  observations  show  that  the  stratified  wake  grows  in  the  hor¬ 
izontal  at  the  exact  same  x 1^3  rate  as  its  unstratified  counterpart  [50,51]. 
Least-square  fits  to  the  SMPM  data  for  both  Fr  =  4  cases  indicate  power  law 
growth  rates  within  the  experimental  uncertainty.  The  corresponding  least 
square  fits  for  the  Fr  =  oo  data  show  power  law  exponents  very  close  to  a 
value  of  1/3.  Throughout  its  entire  evolution,  the  unstratified  wake  maintains 
comparable  values  of  Ly  and  Lz  (not  shown)  which  is  indicative  of  an  ax- 
isymmetric  growth,  as  one  would  expect  [51]. 

Finally,  timeseries  of  the  wake  vertical  lengthscale  Lz  for  the  Fr  =  4  runs 
are  shown  in  figure  8b.  Stratified  wake  experiments  exhibit  a  phase  of  con¬ 
stant  Lz  corresponding  to  the  NEQ  regime  followed  by  a  transition  to  a  f0'6 
growth  rate  in  the  Q2D  regime  [49] .  Both  trends  are  replicated  by  the  SMPM 
data.  The  constant  Lz  phase  does  end  earlier  in  the  Re  =  5  x  103  run.  Nev¬ 
ertheless,  both  Fr  =  4  simulations  transition  into  a  growth  phase  with  power 
law  exponents  that  appear  to  be  closer  to  the  experimental  values  than  those 
of  the  mixed  model  LES  data.  The  slight  discrepancies  in  duration  of  the 
NEQ  regime  and  the  appearance  of  the  breakpoint  when  the  Q2D  growth 
rate  is  established,  should  not  be  grounds  for  questioning  the  reliability  of 
the  spectral  multidomain  model.  As  discussed  in  detail  by  Speckling  [47,49], 
both  the  duration  of  the  NEQ  regime  and  the  transition  into  the  Q2D  phase 
appear  to  be  dependent  on  the  initial  condition  and  specific  forcing  of  the 
flow.  Given  that  (see  §4.2)  the  initial  condition  of  a  stratified  wake  simulation 
at  x/D  =  6  does  not  match  exactly  its  experimental  counterpart  (asymme¬ 
try  in  mean/fluctuating  profiles  zero  density  perturbations  etc.)  these  weak 
discrepancies  are  justifiable. 
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5  Summary  and  Concluding  Remarks 


A  spectral  multidomain  penalty  method  model  has  been  developed  for  the  so¬ 
lution  of  the  three-dimensional  incompressible  Navier-Stokes  equations  under 
the  Boussinesq  approximation  at  high  Re.  A  high  accuracy  splitting  scheme 
is  used  based  on  a  third  order  stiffly  stable  scheme  for  the  non-linear  term  ap¬ 
proximation,  third  order  backward  differentiation  for  the  temporal  derivatives 
and  a  high-order  numerical  boundary  condition  for  the  pressure.  High  spa¬ 
tial  accuracy  is  established  by  Fourier  spectral  discretization  in  the  periodic 
horizontal  and  variable  height  Legendre  spectral  multidomain  discretization 
in  the  vertical,  the  latter  allowing  maximum  flexibility  in  choice  of  boundary 
conditions  and  internal  flow  resolution.  Simulations  of  high  Re  at  resolutions 
of  affordable  cost  are  inherently  under-resolved.  To  overcome  problems  of  nu¬ 
merical  stability  due  to  under-resolution,  moderate  spectral  filtering  is  applied 
in  Fourier  and  Legendre  space  and  the  vertical  discretization  is  supplemented 
with  a  multidomain  penalty  scheme  and  strong  adaptive  interfacial  averag¬ 
ing.  This  is  the  first  time  a  penalty  method  is  applied  to  the  simulation  of 
a  high  Re  incompressible  flow,  with  a  particular  emphasis  on  the  subdomain 
interfaces.  The  penalty  method  is  applied  separately  at  two  levels  in  the  nu¬ 
merical  methodology:  the  non-linear  term  advancement  and  the  viscous  term 
treatment.  Thus,  examination  of  the  internal  high  Re  dynamics  is  possible 
without  full  resolution  of  the  thin  numerical/physical  boundary  layers. 

The  laboratory  flow  chosen  to  validate  the  numerical  model  is  a  momentum¬ 
less  stratified  turbulent  wake,  with  and  without  ambient  stratification,  at  two 
different  values  of  Re.  The  multidomain  results  capture  correctly  the  various 
regimes  of  evolution  and  associated  transition  points  observed  in  stratified 
laboratory  wakes.  The  structure  of  the  vorticity  and  internal  wave  fields  of 
the  stratified  wake,  at  late  and  early-times,  respectively  are  also  well  repro¬ 
duced  by  the  model.  The  unstratified  simulations  exhibit  power  law  exponents 
which  agree  with  similarity  theory  while  the  associated  vortical  structure  is 
comparable  to  that  observed  in  previous  DNS  and  LES.  All  results  indicate 
that  the  model-resolved  three-dimensional  flow  fields  contain  all  the  salient 
wake  dynamics. 

The  spectral  multidomain  model  appears  well-suited  for  the  investigation  of 
localized  stratified  turbulent  flows  at  higher  Re.  The  resolution  used,  if  com¬ 
pared  to  that  required  by  DNS  at  similar  Re,  is  a  factor  of  four  less  in  each 
direction.  Therefore,  high  run  repeatibility  is  possible.  In  regards  to  stratified 
turbulent  wakes,  the  universality  of  scaling  laws  may  be  examined  over  a  broad 
range  of  Re  and  Fr.  The  physical  mechanisms  behind  the  associated  power 
law  exponents  may  be  elucidated  by  investigating  the  interaction  between  vor¬ 
tical  and  internal  wave  modes  in  the  non-equilibrium  (NEQ)  regime  and  the 
late-time  vertical  growth  of  the  wake.  The  flexibility  in  local  flow  resolution  of 
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the  multidomain  technique  motivate  application  of  the  model  to  the  study  of 
the  interaction  of  various  localized  stratified  turbulent  flows  interacting  with 
the  vertical  boundaries.  Examples  include  the  subsurface  currents  generated 
by  internal  wave  radiation  from  the  turbulent  wake,  (a  topic  little  explored  in 
the  existing  literature  [2]),  a  stratified  turbulent  patch  [45]  and  resuspension 
effects  induced  at  the  boundary  layer  by  travelling  solitary  internal  waves  [55] . 
Finally,  the  facility  of  implementation  of  complex  boundary  operators  in  the 
multidomain  penalty  scheme  make  it  amenable  to  the  study  of  geophysical 
flows  with  complex  surface  forcing  (wind  stress/buoyancy  flux). 


6  Acknowledgements 


We  thank  Geoff  Spedding  for  encouragement,  physical  insight  and  reading 
through  an  early  draft  of  this  paper.  We  also  thank  John  Boyd,  Mohammed 
Iskandarani,  Adam  Fincham,  Patrice  Meunier,  Jim  Rottman,  Kraig  Winters 
and  Frank  Tse  for  useful  discussion.  PJD  is  indebted  to  Prof.  John  Diarnes- 
sis,  for  unfaltering  moral  support  and  illuminating  instruction  on  orthogonal 
polynomials.  Jack  Gibbons,  John  Burdette  and  Carl  Gibson  were  instrumen¬ 
tal  in  the  funding  of  this  project.  This  work  was  funded  by  ONR  contract 
N00014-001-0756  administered  by  Dr.  Ron  Joslin.  The  work  of  the  last  au¬ 
thor  (JSH)  was  partially  supported  by  NSF  through  an  NSF  Career  Award 
from  the  Division  of  Mathematical  Sciences,  and  by  the  Alfred  P.  Sloan  Foun¬ 
dation  as  a  Sloan  Research  Fellow.  Computation  for  the  work  described  in 
this  paper  was  supported  by  the  University  of  Southern  California  Center  for 
High  Performance  Computing  and  Communications  (www.usc.edu/hpcc). 


References 

[1]  H.  M.  Blackburn  and  S.  Schmidt.  Spectral  element  filtering  techniques  for  large 
eddy  simulation  with  dynamic  estimation.  J.  Comp.  Phys.,  186:610-629,  2003. 

[2]  V.  G.  Bondur  and  N.  N.  Filatov.  Study  of  physical  processes  in  the  coastal  zone 
for  detecting  anthropogenic  impact  by  means  of  remote  sensing.  In  Proceedings 
of  the  7th  Workshop  on  Physical  processes  in  Natural  Waters ,  pages  98- 
103,  Petrozavodsk,  Russia,  2003. 

[3]  J.  P.  Boyd.  Two  comments  on  filtering  (artificial  viscosity)  for  Chebyshev 
and  Legendre  spectral  and  spectral  element  methods:  Preserving  boundary 
conditions  and  interpretation  of  the  filter  as  a  diffusion.  J.  Comp.  Phys., 
143:283-288,  1998. 

[4]  J.  P.  Boyd.  Chebyshev  and  Fourier  Spectral  Methods.  Dover,  Mineola,  New 
York,  2001. 


25 


[5]  D.  L.  Boyer,  D.  B.  Haidvogel,  and  N.  Perenne.  Laboratory-numerical  model 
comparisons  of  canyon  flows:  A  parameter  study.  J.  Phys.  Oceanogr., 
submitted,  2003. 

[6]  C.  Canuto,  M.  Y.  Hussaini,  A.  Quarteroni,  and  T.  A.  Zang.  Spectral  Methods 
in  Fluid  Dynamics.  Springer- Verlag,  1988. 

[7]  A.  B.  Cortesi,  G.  Yadigaroglu,  and  S.  Banerjee.  Numerical  investigation  of  the 
entrainment  and  mixing  processes  in  neutral  and  stably-stratified  mixing  layers. 
Phys.  Fluids ,  11:162-185,  1999. 

[8]  B.  Costa  and  W.  S.  Don.  On  the  computation  of  high  order  pseudospectral 
derivatives.  Appl.  Num.  Math.,  33:151-159,  2000. 

[9]  M.  O.  Deville,  P.  F.  Fischer,  and  E.  H.  Mund.  High  Order  Methods  for 
Incompressible  Fluid  Flow.  Cambridge  University  Press,  2002. 

[10]  P.  J.  Diamessis  and  K.  K.  Nomura.  The  structure  and  dynamics  of  overturns  in 
stably  stratified  homogeneous  turbulence.  J.  Fluid  Mech.,  499:197-229,  2004. 

[11]  J.  A.  Domaradzki.  An  analytic  Green  function’s  method  in  pseudo-spectral 
Navier-Stokes  solvers  for  boundary  layer  and  channel  flows.  J.  Comp.  Phys., 
88(l):232-242,  1990. 

[12]  D.  G.  Donnnermuth,  J.  W.  Rottman,  G.  E.  Innis,  and  E.  A.  Novikov.  Numerical 
simulation  of  the  wake  of  a  towed  sphere  in  a  weakly  stratified  fluid.  J.  Fluid 
Mech.,  473:83-101,  2002. 

[13]  W.  S.  Don,  D.  Gottlieb,  and  J.  H.  Jung.  A  multidomain  spectral  method 
for  supersonic  reactive  flows.  Technical  Report  BrownSC  2002-06,  Applied 
Mathematics  Dept.,  Brown  University,  2002. 

[14]  W.  S.  Don  and  A.  Solomonoff.  Accuracy  and  speed  in  computing  the  Chebyshev 
collocation  derivative.  SIAM  J.  Sci.  Comput.,  16(6):1253-1268,  1995. 

[15]  D.  Etling.  Turbulence  Collapse  in  Stably  Stratified  Flows:  Application  to 
the  Atmosphere,  pages  1-21.  Clarendon  Press,  Oxford,  1993. 

[16]  A.  Fincham,  T.  Maxworthy,  and  G.  Spedding.  Energy  dissipation  and  vortex 
structure  in  freely  decaying,  stratified  grid  turbulence.  Dyn.  Atmos.  Oceans, 
23:155-169,  1996. 

[17]  P.  F.  Fischer  and  J.  S.  Mullen.  Filtering  techniques  for  complex  geometry  flows. 
Commun.  Numer.  Meth.  Engng.,  15:9-18,  1999. 

[18]  D.  Funaro.  A  multidomain  spectral  approximation  of  elliptic  equations. 
Numerical  Methods  for  Partial  Differential  Equations,  2:187-205,  1986. 

[19]  D.  Funaro.  Domain  decomposition  methods  for  pseudo  spectral  approximations. 
Part  I.  Second  order  equations  in  one  dimension.  Numer.  Math.,  52:329-344, 
1988. 


26 


[20]  R.  P.  Garg,  J.  H.  Ferziger,  S.  G.  Monismith,  and  J.  R.  Koseff.  Stably  stratified 
turbulent  channel  flows.  I.  Stratification  regimes  and  turbulence  suppression 
mechanism.  Phys.  Fluids ,  12:2569-94,  2000. 

[21]  F.  X.  Giraldo,  J.  S.  Hesthaven,  and  T.  Warburton.  Nodal  high-order 
discontinuous  Galerkin  methods  for  the  spherical  shallow  water  equations.  J. 
Comp.  Phys.,  181:499-525,  2002. 

[22]  D.  Gottlieb  and  J.  S.  Hesthaven.  Spectral  methods  for  hyperbolic  problems. 
Journal  of  Computational  and  Applied  Mathematics,  128:83-131,  2001. 

[23]  M.  G.  Gourlay,  S.  C.  Arendt,  D.  C.  Fritts,  and  J.  Werne.  Numerical  modeling 
of  initially  turbulent  wakes  with  net  momentum.  Phys.  Fluids,  13:3783-3802, 
2001. 

[24]  J.  L.  Guermond  and  J.  Shen.  Velocity-correction  projection  methods  for 
incompressible  flows.  SIAM  J.  Numer.  Anal.,  41(1):  112-134,  2003. 

[25]  E.  Hairer,  S.  P.  Nprsett,  and  G.  Wanner.  Solving  Ordinary  Differential 
Equations  I.  Non-stiff  Problems.  Springer  Verlag,  Berlin,  1987. 

[26]  H.  Hanazaki.  A  numerical  study  of  three-dimensional  stratified  flow  past  a 
sphere.  J.  Fluid  Mech.,  192:393-419,  1988. 

[27]  J.  S.  Hesthaven.  A  stable  penalty  method  for  the  compressible  Navier-Stokes 
equations:  I.  Open  boundary  conditions.  SIAM  J.  Sci.  Comput.,  17(3):579- 
612,  1996. 

[28]  J.  S.  Hesthaven.  A  stable  penalty  method  for  the  compressible  Navier-Stokes 
equations:  II.  One-dimensional  domain  decomposition  schemes.  SIAM  J.  Sci. 
Comput.,  18(3):658-685,  1997. 

[29]  J.  S.  Hesthaven.  A  stable  penalty  method  for  the  compressible  Navier-Stokes 
equations:  III.  Multidimensional  domain  decomposition  schemes.  SIAM  J.  Sci. 
Comput.,  20(l):62-93,  1998. 

[30]  S.  E.  Holt,  J.  R.  Koseff,  and  J.  H.  Ferziger.  A  numerical  study  of  the  evolution 
and  structure  of  homogeneous  stably  stratified  sheared  turbulence.  J.  Fluid 
Mech.,  237:499-539,  1992. 

[31]  K.  Horiuti.  A  classification  method  for  vortex  sheet  and  tube  structures  in 
turbulent  flows.  Phys.  Fluids,  13(12):3756-3774,  2001. 

[32]  J.  Imberger.  Transport  processes  in  lakes:  review.  In  R.  Margalef,  editor, 
Limnology  Now:  A  Paradigm  of  Planetary  Problems,  pages  99-193,  North 
Holland,  1994.  Elsevier  Science  Publishers  B.V. 

[33]  M.  Iskandarani,  D.  B.  Haidvogel,  and  J.  P.  Boyd.  A  staggered  spectral  element 
model  with  application  to  the  oceanic  shalow  water  equations.  International 
Journal  for  Numerical  Methods  in  Fluids,  20:393-414,  1995. 

[34]  M.  Iskandarani,  D.  B.  Haidvogel,  and  J.  C.  Levin.  A  three-dimensional  spectral 
element  model  for  the  solution  of  the  hydrostatic  primitive  equations.  J.  Comp. 
Phys.,  186:397-42,  2003. 


27 


[35]  F.G.  Jacobitz,  S.  Sarkar,  and  C.W.  Van  Atta.  Direct  numerical  simulations  of 
the  turbulence  evolution  in  a  uniformly  sheared  and  stably  stratified  flow.  J. 
Fluid  Mech.,  342:231-261,  1997. 

[36]  L.  H.  Kantha  and  C.  A.  Clayson.  Numerical  Models  of  Oceans  and  Oceanic 
Processes.  Academic  Press,  San  Diego,  2000. 

[37]  G.-S.  Karamanos  and  G.  E.  Karniadakis.  A  spectral  vanishing  viscosity  method 
for  large-eddy  simulations.  J.  Comput.  Phys.,  163:22-50,  2000. 

[38]  G.  E.  Karniadakis,  M.  Israeli,  and  S.  A.  Orszag.  High-order  splitting  methods 
for  the  incompressible  Navier-Stokes  equations.  J.  Comp.  Phys.,  97:414-443, 
1991. 

[39]  R.  M.  Kirby  and  G.  E.  Karniadakis.  Coarse  resolution  turbulence  simulation 
with  spectral  vanishing  viscosity-large  eddy  simulations  (SVV-LES).  J.  Fluids 
Eng.,  124:886-891,  2002. 

[40]  P.  K.  Kundu.  Fluid  Mechanics.  Academic  Press,  San  Diego,  2002. 

[41]  J.  G.  Levin,  M.  Iskandarani,  and  D.  B.  Haidvogel.  A  spectral  filtering  procedure 
for  eddy-resolving  simulations  with  a  spectral  element  ocean  model.  J.  Comp. 
Phys.,  137:130-154,  1997. 

[42]  J.  T.  Lin  and  Y.  H.  Pao.  Wakes  in  stratified  fluids.  Ann.  Rev.  Fluid  Mech., 
11:317-338,  1979. 

[43]  J.  S.  Mullen  and  P.  F.  Fischer.  Filter-based  stabilization  of  spectral  element 
methods.  C.  R.  Acad.  Sci.  Paris,  Serie  I  -  Analyse  Numerique,  332:265- 
270,  1997. 

[44]  J.  J.  Riley  and  M.  P.  Lelong.  Fluid  motion  in  the  presence  of  strong 
stratification.  Annual  Review  of  Fluid  Mechanics,  32:613-657,  2000. 

[45]  J.  J.  Riley  and  R.  W.  Metcalfe.  Direct  numerical  simulations  of  turbulent 
patches  in  stably-stratified  fluids.  In  Stratified  Flows  I,  California  Institute  of 
Technology,  Pasadena,  1987.  Proc.  3rd  Int.  Symp.  on  Stratified  Flows. 

[46]  D.  N.  Slinn  and  J.  J  Riley.  A  model  for  the  simulation  of  turbulent  boundary 
layers  in  an  incompressible  stratified  flow.  J.  Comp.  Phys.,  144:550-602,  1998. 

[47]  G.  R.  Spedding.  The  evolution  of  initially  turbulent  bluff-body  wakes  at  high 
internal  Froude  number.  J.  Fluid  Mech.,  337:283-301,  1997. 

[48]  G.  R.  Spedding.  Anisotropy  in  turbulence  profiles  of  stratified  wakes.  Phys. 
Fluids,  13(8):2361-2372,  2001. 

[49]  G.  R.  Spedding.  Vertical  structure  in  stratified  wakes  with  high  initial  Froude 
number.  J.  Fluid  Mech.,  454:71-112,  2002. 

[50]  G.  R.  Spedding,  F.  K.  Browand,  and  A.  M.  Fincham.  Turbulence,  similarity 
scaling  and  vortex  geometry  in  the  wake  of  a  towed  sphere  in  a  stably  stratified 
fluid.  J.  Fluid  Mech.,  314:53-103,  1996. 


28 


[51]  H.  Tennekes  and  J.  L.  Lumley.  A  first  course  in  turbulence.  The  MIT  press, 
1972. 

[52]  S.  A.  Thorpe.  Turbulence  in  the  stratified  and  rotating  world  ocean.  Theor. 
Comput.  Fluid  Dyn.,  11:171-181,  1998. 

[53]  J.  Trujillo  and  G.  E.  Karniadakis.  A  penalty  method  for  the  vorticity-velocity 
formulation.  J.  Comp.  Phys.,  149:32-58,  1999. 

[54]  K.  L.  Tse,  A.  Mahalov,  B.  Nicolaenko,  and  H.  J.  S.  Fernando.  Quasi-equilibrium 
dynamics  of  shear-stratified  turbulence  in  a  model  tropospheric  jet.  J.  Fluid 
Mech.,  496:73-103,  2003. 

[55]  B.  J.  Wang  and  L.  G.  Redekopp.  Long  internal  waves  in  shear  flows: 
Topographic  resonance  and  wave-induced  global  instability.  Dyn.  Atmos. 
Oceans ,  33:263-302,  2001. 

[56]  K.  B.  Winters,  J.  McKinnon,  and  B.  Mills.  A  spectral  model  for  process  studies 
of  density  stratified  flows.  J.  of  Atmos.  Ocean.  Techn.  (submitted),  2003. 


29 


u 


Fig.  1.  Computational  domain  for  the  simulation  of  a  mid-to-late  time  momentlum- 
less  stratified  turbulent  wake.  The  wake  was  originally  generated  by  a  sphere  of 
diameter  D,  towed  with  a  velocity  U,  which  however  is  not  present  in  the  compu¬ 
tational  domain.  The  domain  dimensions  are  Lx  x  Ly  x  Lz.  Consistent  with 
the  salt-stratified  water  tank  the  domain  has  a  solid  wall  bottom  and  free-slip  top 
and  when  needed,  a  linear  stable  density  stratification  may  be  imposed.  Shown  are 
also  the  horizontal  and  vertical  centerplane  where  the  experimental  data  used  in 
the  initial  condition  are  sampled. 
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Fig.  2.  Oxz  section  of  the  numerical  grid  employed  in  this  study.  The  horizontal 
direction  employs  a  spectral  Fourier  discretization  and  uniform  grid.  In  the  vertical 
a  Legendre  spectral  multidomain  discretization  is  used  with  M  =  5  subdomains 
of  order  of  approximation  N  and  a  local  Gauss-Lobatto-Legendre  grid.  Shown  is 
the  vertical  grid  point  distribution  for  Re  =  5  x  103  where  N  =  32.  Re  =  2  x  104 
employs  N  =  48.  The  dashed  horizontal  lines  on  the  right  of  the  figure  indicate  the 
position  of  the  subdomain  interfaces.  Subdomain  origins  are  at  z/D  =  0,  3.2,  5.2, 6.8 
and  8.8. 
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Fig.  3.  Normalized  one-dimensional  spectra  of  turbulent  kinetic  energy  calculated 
for  Re  =  5  x  103  and  2  x  10.  The  streamwise  spectra  are  averages  of  estimates  taken 
at  all  spanwise  locations  (including  the  wake  far  field)  on  the  horizontal  centerplane. 
Spectra  are  shown  at  initialization  of  the  actual  run  as  well  as  intermediate  and  late 
times. 
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(a) 


(b) 

Fig.  4.  Isosurfaces  of  vertical  vorticity,  ujz  in  a  y  and  z-truncated  domain  for  the 
Fr  =  4  simulations.  Also  shown  are  contours  of  wz  at  horizontal  centerplane 
(z  =  Lz/2).  (a)  Re  =  5  x  103  ,  (b)  Re  =  2  x  104.  Isosurface  values  are  \l2\u:z\max 
and  —^/2\ioz\max  (light  and  dark  color,  respectively).  Minimum/  maximum  contour 
values  on  centerplane  slice  are  taken  as  half  the  minimum/  maximum  uz  value  in 
the  flow. 
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Fig.  5.  Isosurfaces  of  vorticity  magnitude  |u;|  in  a  y  and  z-truncated  domain  for  the 
Fr  =  oo  simulations,  (a)  Re  =  5  x  103  ,  (b)  Re  =  2  x  104.  Isosurface  value  is 
1/3M  max • 


(a)  (b) 


Fig.  6.  Contours  of  horizontal  velocity  divergence,  A~  at  the  vertical  centerplane 
of  the  flow  ( x  =  0.5 Lx),  for  Fr  =  4  sampled  at  time  Nt  =  15.  (a)  Re  =  5  •  103, 
(b)  Re  =  2  •  104.  Minimum/maximum  contour  values  are  taken  as  half  the  mini¬ 
mum/maximum  value  of  A,  in  the  flow.  The  arrows  indicate  direction  of  internal 
wave  propagation  inclined  at  45°  with  respect  to  the  vertical. 
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Fig.  7.  Timeseries  of  the  appropriate  scaled  maximum  centerline  mean  velocity  Uq. 
Comparison  of  Fr  =  4  and  oo  SMPM  results  with  data  from  experiments,  DNS  and 
mixed  model  LES.  The  Nt  and  ( U/Uq )  •  Fr 2/3  scaling  for  the  abscissa  and  ordi¬ 
nate,  respectively,  is  valid  for  all  the  stratified  data.  The  x/ (2D)  and  (U /Uq)  ■  Fr 4/3 
scaling  is  valid  for  comparing  Fr  =  4  to  Fr  =  oo  results.  The  dotted  lines  are 
the  laboratory  data  of  Spedding  [50,47]  where  the  middle  line  represents  the  mean 
value  and  the  top/bottom  lines  correspond  to  the  maximum  deviation  in  the  exper¬ 
imentally  observed  values.  The  vertical  lines  delineate  the  three  distinct  regimes  of 
wake  evolution  [47].  Also  shown  is  the  x~2^  power- law  predicted  by  self-similarity 
theory. 
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Fig.  8.  For  caption  see  following  page. 
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Fig.  8.  Timeseries  of  appropriately  scaled  mean  velocity  profile  lengthscales.  The 
respective  SMPM  results  are  compared  to  the  relevant  Spedding  laboratory  data, 
DNS  and  mixed  model  LES.  Dotted  lines  represent  corresponding  Spedding  exper¬ 
imental  data  [50,47,49]  as  explained  in  figure  7.  (a)  SMPM  results  for  horizontal 
wake  lengthscale,  Ly  at  both  Fr  =  4  and  oo.  The  Fr  =  4  Spedding  experimental 
data  evolves  with  the  x1/3  power  law  predicted  by  similarity  theory.  The  horizontal 
variable  is  x/D.  Lines  with  full  triangles  are  normalized  Fr  =  oo  values  of  vertical 
lengthscale  Lz ■  (b)  SMPM  results  for  vertical  wake  lengthscale,  Lz  at  Fr  =  4.  The 
horizontal  variable  is  Nt. 
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